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ABSTRACT 

We present an interface between the (magneto-) hydrodynamics code PLUTO and the plasma simulation and spectral synthesis code 
CLOUDY. By combining these codes, we constructed a new photoionization hydrodynamics solver: The PLUTO-CLOUDY Interface 
(TPCI), which is well suited to simulate photoevaporative flows under strong irradiation. The code includes the electromagnetic 
spectrum from X-rays to the radio range and solves the photoionization and chemical network of the 30 lightest elements. TPCI follows 
an iterative numerical scheme: First, the equilibrium state of the medium is solved for a given radiation field hy CLOUDY, resulting 
in a net radiative heating or cooling. In the second step, the latter influences the (magneto-) hydrodynamic evolution calculated hy 
PLUTO. Here, we validated the one-dimensional version of the code on the basis of four test problems: Photoevaporation of a cool 
hydrogen cloud, cooling of coronal plasma, formation of a Strdmgren sphere, and the evaporating atmosphere of a hot Jupiter. This 
combination of an equilibrium photoionization solver with a general MHD code provides an advanced simulation tool applicable to a 
variety of astrophysical problems. 
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1. Introduction 


Hydrodynamic flows powered by strong irradiation (photoevap¬ 
oration) can be found throughout the Universe from the evap¬ 
oration of cosmological minihaloes ( [Shapiro et al.||20()4|l a nd 
the formation and evolution of Hit regions (O’Dell 200 l| l to 
the evaporation of circumstellar disks ( Owen et al.||2()10]r~^e 
discovery of expanding hot-Jupiter atmospheres ( jVidal-Madjar] 
et al.|2007j l has added another environment possibly dominated 
by photoevaporative mass loss ( Lammer et al.|2005| . Although 
photoevaporative phenomena occur on widely different spacial 
scales, the essential physical processes are similar: the absorp¬ 
tion of ionizing radiation causes an increase in temperature and 
pressure, which results in a hydrodynamic wind. 

Coupled photoionization and hydrodynamic simulations are 
essential for the progress of research in these fields. We created 
a new interface between an extensive equilibrium solver for the 
state of a gas or plasma under strong irradiation (CLOUDY) and 
a 3D (magneto-) hydrodynamics simulation code (PLUTO), the 
PLUTO-CLOUDY Inter face (TPCI). The codes and the inter¬ 
face are publicly available^ The interface can be used to study 
photoevaporative processes under a wide range of conditions. 
However, we designed the code to study the environment of hot- 
Jupiter atmospheres, and a short introduction provides an idea of 
our requirements for the new simulation tool. 

The photoevaporation of hot-Jupiter atmospheres is the re¬ 
sult of a complex interplay of various physical processes. The ab¬ 
sorption of extreme ultraviolet radiation (EUV, A = 100-912 A) 


* http://www.nublado.org/ 

^ http://plutocode.ph.unito.it/ 

^ http://www.hs.uni-hamburg.de/DE/Ins/Per/Salz/ 


heats the upper atmosphere of planets and in extreme cases trig¬ 
gers expansion and evaporation ( [Watson et al.|198T]). Our focus 
lies on active host stars such as Corot-2 ( [Schroter et al.||2011| l, 
where the effect of X-rays (T < 100 A) on the mass-loss rate can¬ 
not be neglected ( Cecchi-Pestellini et al.[[2006 i. With standard 
abundances, metals dominate the absorption of X-rays (Morri 


& McCammon||1983[ l and are in many cases important for 
the cooling of gases that are devoid of molecules and dust (e.g., 
in Hit regions, see Osterbrock & Ferland|[2006) l. Thermal con¬ 
duction can be a major heat source if strong temperature gra¬ 
dients occur in the atmospheres ( [Watson et al.lllOST) . Further¬ 
more, hot Jupiters are most likely tidally locked, and eventually 
multidimensional simulations are needed to solve the disparate 
flows from the day and night side of the planets. Finally, plan¬ 
etary magnetic fields ([Trammell et al.||2011[) or the interaction 
with the stellar wind ( [Tremblin & Chiang||2013[ l also affect the 
evaporation process. 


In this paper, we describe and validate the one-dimensional 
version of our numerical scheme without focusing on a particular 
science application. We start by describing our new coupled pho¬ 
toionization hydrodynamics simulation scheme (Sect. HI} and 
comparing it to available numerical codes (Sect. |2.2[ l. The the¬ 
oretical background and additional aspects of the interface are 
explained in Sects. [Z^ to |2.8[ To verify the results of the code 
and demonstrate the range of applications, four problems were 
drawn from different fields and were solved with TPCI (Sect.j^. 
Our last test case is a simplified simulation of the escaping at¬ 
mosphere of the hot Jupiter HD 209458 b, which we compare to 
previous studies of the same system (Sect. |3.4[ l. In Sect. we 
discuss the results and indicate future applications for TPCI. 
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Fig. 1. Flowchart of the PLUTO-CLOUDY Interface. Affiliation of modules is indicated by colored boxes. PLUTO solves the hydrodynamic 
evolution and CLOUDY the ionization and chemical equilibrium state of the gas under strong irradiation. Communication is achieved through the 
interface, which also handles the interpolation between the two different grid structures. 


2. The PLUTO-CLOUDY Interface - TPCI 

2.1. Overview 

We introduce an interface between two popular codes: the 
(magneto-) hydrodynamics code PLUTO ( [Mignone et al.|20T2 
and the photoionization simulation code CLOUDY ( |Ferland 
et al. 2013|l. The interface allows studying steady-state and 


slowly evolving photoevaporative flows, in which the medium 
is in ionization and chemical equilibrium at all times during the 
simulation. 

The numerical scheme follows a clear division between the 
solution of the (magneto-) hydrodynamic equations in PLUTO 
and the solution of the ionization and chemical equilibrium in 
CLOUDY (see Fig. [^I. Basically, the photoionization solver is 
called by the hydrodynamic code with the density and temper¬ 
ature distributions and solves the microphysical state of the gas 
under the given conditions and with the specified irradiation. It 
returns the radiative heating or cooling distribution, which then 
affects the hydrodynamic evolution. The two steps are alternated 
during the simulation. 

PLUTO is a freely distributed modular code for solv¬ 
ing HD/MHD equations ( [Mignone et al.|[2007| 2012| l. It is a 
Godunov-type code that computes intercellular fluxes by solv¬ 
ing the Riemann problem at cell interfaces. Parabolic terms in 
the differential equations such as viscosity or conductivity can be 
included by operator splitting, which is also the case for source 
terms such as optically thin cooling or gravity. For TPCI we used 
the static grid version of the code. 

CLOUDY is a photoionization and microphysics equilibrium 
solver for static structures, irradiated by an arbitrary source ( |Fer- 
jland et al.|1998| 2013| l. It accounts for the electromagnetic spec¬ 
trum from hard X-rays to the radio regime. The 30 lightest ele¬ 
ments, from hydrogen to zinc, are included, and the solver bal¬ 
ances radiative and collisional ionization and recombination, but 
also more complicated physical processes such as inner shell 
ionization or charge exchang^ In TPCI CLOUDY is called with 
a fixed temperature and density structure and solves the ioniza¬ 
tion and chemical equilibrium for the given conditions, which 
results in an imbalance of radiative heating and cooling. The pro¬ 
gram has been thoroughly tested to approach LTE at high den¬ 
sities and can be used within the temperature range from 3 K to 
10'° K and for number densities of up to 10^^ cm“^. In cold re¬ 
gions CLOUDY approaches the fully molecular limit, including 


We sometimes us the generic term microphysics, meaning all pro¬ 
cesses that affect the equilibrium state in CLOUDY. For a full list we 
refer to the CLOUDY publications. 
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about 10^ chemical reactions mostly from the UMIST database 
( |Le Teuff et al.||2000l l. The ionization and chemical network is 
solved completely within CLOUDY, and only the mean molecu¬ 
lar weight of the medium is passed to PLUTO. Since CLOUDY 
is a static equilibrium solver, advection of species is in general 
not included in TPCI. However, in ID simulations with with 
a bulk flow antiparallel to the direction of irradiation, advec- 
tive effects on the ionization equilibrium can be included via a 
CLOUDY internal iterative solver. Diffusion of elements is not 
included in this approach. 

The coupled simulation codes use independent grid struc¬ 
tures, and communication is achieved by linear interpolation 
on tabulated structures. TPCI is capable of performing true 3D 
MHD simulations coupled to pseudo-3D photoionization simu¬ 
lations along ID parallel slices through the computational do¬ 
main. One-dimensional simulations with the interface are serial, 
while 2D or 3D simulations are parallelized through the message 
passing interface (MPI). Here, we restrict the verification to ID 
problems without magnetic fields. However, multidimensional 
simulations are not treated differently in the interface. 

In the design of the interface we followed the paradigm to ap¬ 
ply only minimal changes to both programs so as not to interfere 
with the well-tested numerical schemes while exploiting most of 
the capabilities of both codes. To this end, CLOUDY was used 
as an external library for solving the equilibrium state of the gas. 
Several input and output extensions were necessary to enable 
the use of CLOUDY within TPCI: we adapted the code to ac¬ 
cept tabulated density, temperature, and velocity structures, and 
the methods for retrieving the results from CLOUDY were ex¬ 
panded to include all necessary data in the communication with 
PLUTO. 

The setup of TPCI follows the setup of the individual codes. 
While the usual input script for CLOUDY simulations is speci¬ 
fied directly in the interface, the input files for the PLUTO code 
have not been changed. For the output of the hydrodynamic vari¬ 
ables the standard PLUTO output options are available. Through 
the CLOUDY output much additional information is available, 
such as the number densities of all ions and molecules, the im¬ 
portant radiative heating and cooling agents with contribution to 
the total heating or cooling rate, and the irradiating, transmitted, 
and reflected spectra. In general, anybody who is familiar with 
the two simulation programs individually will be able to use the 
interface with only little additional introduction. 
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2.2. Other photoionization hydrodynamics simuiation codes 

For reference, we provide a short overview of some available 
radiation hydrodynamics codes and compare individual aspects 
of the numerical schemes with our new interface. 

For instance, Owen et^(2010) coupled the 3D photoion¬ 


ization and radiative transfer code MOCASSIN (Ercolano et al. 


2003 I with the hydrodynamic code ZEUS-2D ([Stone & Norman 


1992 I. MOCASSIN and CLOUDY are similarly extensive equi¬ 


librium photoionization solvers, but the temperature parameteri¬ 
zation used by [Owen et al.j is only valid for X-ray heating. 

Another example for coupling a microphysical equilibrium 
solver to a hydrodynamic simulation is the ionization module 
for the ELASH code ( jEryxell et al.|200()|) presented byjRijkhorstj 
et al.| ( 2006 1 and further improved by Peters et al. (20101. The 3D 


radiative transfer method is highly efficient in simulations with 
adaptive mesh refinement on distributed systems, but is com¬ 
putationally more demanding than the pseudo-3D scheme used 
here. A similarly advanced parallel radiative transfer method was 
introduced by [Wise & Abelj (|2011|l into the ENZO code (jBryanj 
|& Norman||1997||0’Shea et al.||2004j l; it is called MORRAY. 
The non-equilibrium chemistry solver is restricted to hydrogen 
and helium, however ([The Enzo Collaboration et al.|2013]l. 


Shapiro et al.j (20041 and predating publications have ex¬ 
tended the hydrodynamics code CORAL (jRaga et al. |1995|l 


to include radiative transfer and non-equilibrium photoioniza¬ 
tion of hydrogen, helium, and metals. The scheme uses a simi¬ 
lar pseudo-multidimensional radiative transfer method, but ne¬ 
glects X-rays, which is one of our main interests. X-rays are 
also mostly neglected in the numerical schemes, which have 
been specifically designed to simulate escaping hot-Jupiter at¬ 
mospheres, and the authors focus exclusively on ID simulations 


(e.g., [Yelle 

[20041 [Tian et al.[ 20051 

Garcia Munoz [2007 [ Penz 

[et al.|2008| 

Murray-Clay et al.|2009[[Koskinen et al.|2013[l. 


In comparison, only TPCI solves our need for a photoion¬ 
ization hydrodynamics solver including hydrogen, helium, and 
metals as well as the absorption of EUV and X-ray emission. 


2.3. Hydrodynamic simuiation — PLUTO 


In the following we present only the equations that are most 
relevant for the theoretical background of TPCI. The hydrody¬ 
namic simulation part is solved by the code PLUTO ( [Mignone 
et al.|20d7| 2012[ l. Eor a single fluid with mass density p, pres¬ 
sure p, and velocity v, the hydrodynamic equations in the quasi¬ 
conservative form are given by 


f F = v.n + s 


P ■ 


( 1 ) 


Here U = (p,pv, £) denotes the vector of conserved variables: 
mass, momentum, and total energy density, which is given by 
£ = p/(r - 1) -H l/2pv^, where F (= 5/3) is the specific heat 
ratio of a monatomic gas. F = (pv,pvv + p,{& + p)\) is the hy¬ 
perbolic flux tensor, and H is the parabolic flux tensor, which in 
our case only contains the conductive energy flux Fc. Sp repre¬ 
sents source terms such as gravity Sg = (0,paG,pv ■ Bg) with 
the gravitational acceleration Bg, and the radiative source term 
Sr = (0,0, (Gr - Lr)), where Gr and Lr are the radiative heat¬ 
ing and cooling rates. 

We follow a standard operator split scheme and include the 
radiative heating or cooling rate between the hydrodynamic in¬ 
tegration steps. First, the homogeneous left-hand side of Eq. [T] 
namely the Euler equations, are solved. Second, PLUTO com¬ 
putes the included parabolic and source terms. Third, CLOUDY 


is called to compute the radiative source term (see Sect. \2.4} , 
which subsequently affects the internal energy through the ther¬ 
mal pressure; 

p«+i H-Af"(r-1)5;^, (2) 

where At" denotes the explicit time step at integration step n. 

To minimize the computational effort, the radiative source 
term is not updated at every time step, but only if a change in 
either pressure or density by more than a user defined change 
factor triggers a call to the photoionization solver. A value 
of 1% for the change factor is a good compromise between 
accuracy and convergence speed for an initial settling phase of 
steady-state solutions. The accuracy can be increased by restart¬ 
ing the code with a smaller change factor. 


2.4. Photoionization soiver — CLOUDY 


TPCI uses CLOUDY to solve the photoionization and the equi¬ 
librium state of the medium. We provide a compact overview of 
the most important processes for our applications in CLOUDY; 
for a more general introduction see [Eerland et al. ( [1998 2013[ l. 
The code steps through a one-dimensional slice of the computa¬ 
tional domain with an adaptive step width using linear interpo¬ 
lation on the tabulated density and temperature stractures passed 
by PLUTO. In each cell it solves the local equilibrium state 
by taking into account ionization and recombination processes, 
chemical reactions, and atomic level transitions of all included 
elements. 

CLOUDY separates the radiative transfer of continuum and 
line radiation. The continuum is diminished by absorption, 


/ = /O exp (-dTabs) , 


(3) 


where I is the intensity after the initial intensity Iq passes a slab 
of gas with the total optical depth drabs- The total optical depth 
contains all continuous opacity sources included in CLOUDY, 
but in the EUV range it is dominated by ionization processes of 
neutral hydrogen. The continuum is furthermore diminished by 
scattering ( [Lightman & White|1988[ l: 

/ = fo(l + ^dTscat) ■ (4) 


The total scattering optical depth Tscat also contains the opacities 
of several processes, of which Rayleigh scattering in the Lyman 
line wings dominates in the UV range. 

The escape probability mechanism ( [Castor[ri970[ [Elitzur[ 
1982[ l is used to approximate radiative transfer effects. The lo¬ 
cal mean intensity J, needed to solve the rate equations, is then 
given by 


J^S(l- Peso) . 


(5) 


Here S is the source function and Pgsc the escape probability, 
which only depends on the optical depth. Solving the local equi¬ 
librium state is thereby decoupled from the radiative transfer 
equation. This approximation is exact in a single cell, where 
the source function is constant ( [Elitzur[1992[ l. CLOUDY ensures 
this by choosing the adaptive grid resolution so that neither den¬ 
sity, temperature, heating, or other properties change strongly 
within one cell. However, it is necessary to step through the com¬ 
putational domain at least twice to obtain a reasonable estimate 
for the optical depth in the computation of the escape probabil¬ 
ity. This iteration is needed because radiation can escape in both 
directions, that is, away and toward the source. 
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CLOUDY can include the 30 lightest elements with any 
abundances. These species are only present in the photoioniza¬ 
tion solver, only the mean molecular weight is passed to the hy¬ 
drodynamic single fluid simulation. The ionization ladder of ev¬ 
ery element is computed by balancing the rate equations 

T) - n{X‘) (r(Y') -H nMX‘^\T)) = 0 

( 6 ) 

for each ionization stage. Here n{X‘) is the number density of 
the element X in the ionization stage i, is the electron num¬ 
ber density, a{X\ T) is the temperature-dependent recombina¬ 
tion rate from to X' summed over all levels, and r(Y') is the 
photoionization rate given by 



Here Jy is the mean intensity at frequency v, ay(X‘) is the pho¬ 
toionization cross-section, and vq,, is the threshold frequency 
above which photoionization is possible. For simplicity we do 
not show more complex physical processes such as secondary 
ionizations or charge exchange in Eq. which are also solved 
by CLOUDY, however. 

For the hydrodynamic part of the simulation we need the ra¬ 
diative heating or cooling rates. While the heating rate Gr is a 
sum of many processes, the main heat source is usually given by 
photoionization: 


AjrJ 

I —^h{v-voj)ay(X‘)dv, (8) 

Jvo., 

the sum includes all species. Locally other heat sources such 
as line absorption, charge exchange, and chemical reactions can 
dominate the heating rate. 

Many processes contribute to the cooling rate Lr of a gas: 
recombination processes, free-free emission of electrons in the 
held of ions, line radiation, free-bound transition of H“ , and 
chemical reactions, to name a few. Cooling due to recombina¬ 
tions is given by 

Ll = nMX‘^^)kTI3{r, , (9) 

where yS(Y', T) stands for the effective recombination coefficient, 
averaged over the kinetic energy of the electron gas. The escape 
probability Pgsc ensures that recombinations only contribute to 
the cooling rate if the emitted photon escapes the medium. An¬ 
other main cooling agent is line radiation of collisionaly excited 
elements: 


Gr = 2n(Y') 

X' 


Ll = Me (niClu - MuCul) hv\u ■ 


( 10 ) 


Here ciu denotes the collision rate between the lower and up¬ 
per levels of the atom. The difference of collisional excitation 
and de-excitation is emitted radiatively and leads to cooling. 
The level populations of each ionization stage of every element 
are solved by balancing the equilibrium radiative and collisional 
transition rates in model atoms of differing complexif y. Thus, 


CLOUDY includes a large number of lines (10^ 
let al.|2013| l. 


10 '’ 


Ferland 
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2.5. Time step control and relevant timescales 

In TPCI simulations, restrictions due to hydrodynamic and mi¬ 
crophysical timescales must be observed. PLUTO chooses the 
next hydrodynamic time step based on the Courant-Friedrichs- 
Levy (CFL) condition ( Courant et al.|1928l l, which basically re¬ 
stricts disturbances in the medium to propagate by less than one 
cell width per time step. Radiative heating and cooling addition¬ 
ally limits the hydrodynamic time step in TPCI, meaning that 
the energy loss due to radiative cooling cannot exceed the inter¬ 
nal energy. It is efficient to limit the next time step depending on 
the fraction of the energy content and the radiative source term 
( |Frank & Mellema|1994| l: 


Atl 


TPCI 


cr 


(F-l)5” ’ 


( 11 ) 


and the next hydrodynamic time step is then given by 


= min [AfpLu.pQ, Af^p^ . (12) 

Here cr is the user-defined maximum fractional variation param¬ 
eter. The default value is Cr = 0.1, which restricts the change in 
internal energy due to radiative heating or cooling to be lower 
than 10% in every time step. 

CLOUDY solves the equilibrium state of the medium based 
on a diversity of processes, each with a characteristic timescale, 
and the use of TPCI implies that the microphysical timescales 
are shorter than the hydrodynamic timescale. In gases with a 
temperature of about 10 000 K, hydrogen recombination usually 
is the microphysical process with the longest timescale (|Ferland| 
[T9^ : 

Lrec = ,1. ^ 1.5 X lO*" . (13) 

aA(Te)ne 


TPCI locally checks whether recombination is slower than ad- 
vection and issues a warning. Equation 13 holds only for ion¬ 
ized regions; especially molecular reactions occur on longer 
timescales. The CLOUDY output provides a manual method to 
control the longest timescales and, thus, check the validity of the 
approach. 

In ID simulations, the advection of species can be included 
through a steady-state solver in CLOUDY (see Sect. 2.6 1 . With 
this scheme, steady flows can be simulated in which the advec- 
tive timescales are shorter than microphysical timescales. Nev¬ 
ertheless, the assumption of photoionization and chemical equi¬ 
librium is a strong simplification, and the validity must be con¬ 
sidered for every application. 


2.6. Advection of species 

The impact of advection due to a bulk flow of the medium can 
be included in the simulation by using an iterative steady state 
solver included in CLOUDY ( |Henney et al.||2005| ). The solver 
is available for simulations in which the medium flows toward 
the source of ionizing radiation. The solution of the equilibrium 
state then includes source and sink terms of all species, induced 
by material flows: 

V ■ (n(Y')v) = nv ■ V (n(Y')/n) . (14) 


In this case, PLUTO passes the velocity in addition to the density 
and temperature to the photoionization solver. CLOUDY com¬ 


putes the advective source terms using the relative density (Hen- 
|ney et al.|2005|l, which is an advective scalar. 
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The steady-state solver in CLOUDY includes advection in 
an iterative manner, and the user must control the convergence 
of the solution. The module defines the advection length Az, a 
measure for the resolution of advection processes. It is refined 
if the advective solution has converged for the current advection 
length. The progress is controlled by following the discretization 
error and the convergence error (see |Henney et al.||20d5] for a 
detailed description). The convergence error is given by 


e 


2 

1 


Z 


/ 


\ 


nJ{X‘) - nJ-\X‘) 
Aziv 



/ 


(15) 


where m is the iteration number and j is the zone number; the 
sum extends over all elements and ionization stages as well as 
over all zones. The error approaches zero as the changes in num¬ 
ber densities of consecutive iterations decreases. The discretiza¬ 
tion error is defined as 


= Z 

x‘j 


nJiX') ■ 




n’J(X‘) ■ 




AzIv 


Az/2v 


(16) 


(j-Az) is a symbolic index referencing the density of the species 
at the distance of the advection length. The error describes the 
differences of the local, advective source or sink terms to the 
value if the advection length is refined by a factor of two. The 
sum again runs over all species and zones. 

The discretization error depends on the resolution in the pho¬ 
toionization solver. This resolution can be up to ten times higher 
than in the hydrodynamic part of the simulation. For TPCI it 
is more important to check the advection length than the actual 
value of the discretization error. The advection length is refined 
ifef < 0.1 e|. It should generally be decreased to the hydro- 
dynamic grid spacing, so that advective effects have the same 
accuracy as the rest of the simulation. 

The described procedure does not follow our split scheme to 
solve hydrodynamic effects exclusively in PLUTO. However, at 
the moment we cannot solve the advection of species in PLUTO. 
It is possible to include advective scalars for the relative den¬ 
sity of all species in PLUTO, but then CLOUDY would have 
to be called with the densities of all 30 elements and ioniza¬ 
tion stages, which is not implemented. The steady-state solver in 
CLOUDY is computationally demanding. One advection solu¬ 
tion may need up to 100 iterations to converge because advective 


effects must be swept through the computational domain (Hen- 
ney et al. 2005| l. Thus, including advection slows simulations 
down, and when searching for advective steady-state solutions, 
it is advisable to use the non-advective code for an initial settling 
of the hydrodynamic structure and then restart the simulation in¬ 
cluding advection to let the structure evolve to the final steady 
state. This procedure has been followed in the two test problems, 
which include advection of species (see Sects. [TT|and|3.4|l. 


2.7. Radiative acceleration 

In addition to affecting the temperature of the gas, the absorp¬ 
tion of radiation also exerts a pressure on the medium. CLOUDY 
records the radiative pressure produced by the absorption of ra¬ 
diation from the central source. The local radiative acceleration 
can be directly accessed after a CLOUDY call. This accelera¬ 
tion is handled like any external force by the PLUTO code (see 

Eq.g. 


2.8. Pseudo-multidimensional simulations 


Although here we focus on ID simulations, TPCI can be used 
in two or three dimensions on a Cartesian grid with irradiation 
along the x-axis. The radiative solution is then simply split into 
parallel ID simulations, which is called a pseudo-3D radiative 
transfer. There is no true 3D radiative transfer in CLOUDY, but 
the coupling of 3D radiative transfer to hydrodynamics simu¬ 
lations is in many cases still computationally prohibitive today 
(e.g., |Owen et aH2010| l. In problems with a strong directional¬ 
ity in the radiative field, the pseudo-3D approximation is valid, 
unless there are shadow-casting objects ( Morisset et al.|[2005| l. 
However, in the case of strong shadows, the induced error can 
also be neglected if the energy transfer into the shadow region is 
hydrodynamically dominated. 

The ID version of the interface is a serial code because 
CLOUDY requires the complete density structure to solve the 
plane-parallel radiative transfer through the computational do¬ 
main. PLUTO is, however, capable of 3D parallel computations, 
which is also enabled in TPCI by not decomposing the domain 
along the x-axis. Each thread then deals with the full ID struc¬ 
ture along the x-axis for a given number of y/z-grid-points. Af¬ 
ter the complete domain has been processed, the hydrodynamic 
evolution continues until the next call to CLOUDY is necessary. 
Eor this approach the computational domain should be decom¬ 
posed so that the available processors have the same number of 
y/z-points to minimize idle times. 

This pseudo-3D scheme provides a convenient possibility 
to allow parallel simulations with low complexity. In the best 
case, a multidimensional simulation will have a similar execu¬ 
tion time as a ID simulation if every y/z-grid-point uses an in¬ 
dividual processor. The communication overhead is small, since 
the main computational load usually comes from the indepen¬ 
dent CLOUDY calls. 


3. Verification 

In the following we present TPCI simulations of four test prob¬ 
lems. The first two simulations verify the implementation of our 
numerical scheme, each with a different focus. The third test 
problem is included to clarify the limits of validity of TPCI. The 
fourth test case is a simplified simulation of a hot-Jupiter atmo¬ 
sphere, which already demonstrates the capabilities of TPCI in 
this field. The setup of each simulation is given in Table 


3.1. Weak D-type ionization front 

Eor the initial verification, we investigated the effects of steady- 
state flows in ID simulations of ionization fronts. Such flows 
appear when molecular clouds are evaporated by the ionizing 
emission of newly formed stellar clusters. The gas is heated 
through ionization processes, and the cloud slowly evaporates. 
If the evaporation proceeds as a subsonic wind from a dense re¬ 
gion, it is called a weak D-type ionization front, in contrast to 
strong ionization fronts, which involve transonic flows. (Eor rar¬ 
efied R-type ionization fronts see Sect. 3.3 ) 

The results from our TPCI simulation are compared with an 
independent CLOUDY simulation, using the included steady- 
state solver, which was s pecifically des i gned t o solve these pho- 
toevaporative flows (see Henney et al. (2005 ( ^ This is a pow¬ 
erful test for TPCI because both simulations rely on the same 


^ The CLOUDY internal, iterative, steady-state hydrodynamics solver 
is restricted to specific flows and geometries. TPCI is more versatile and 
solves the dynamic evolution of the studied flows. 
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Table 1. Simulations for the verification of TPCI 


Name 

Numerical setup 

(time stepping, interpolation, 

Riemann solver) 

BC 

left 

pvp 

right 

pvp 

Grid 

geometry 

grid points 

Irradiation 
SED shape 

ioniz. radiation field 
(A < 912 A) 

D-front 

RK3, WEN03, hllc 
-i-advection 

f f o 

O f 0 

Cartesian 

200 uni. 

50000KBB 

(pjj = lO" cm“^ s“' 

c-shock 

RK3, LINEAR, two_shock 

ooo 

OOO 

Cartesian 

400 uni. 

— 

— 

R-front 

RK3, LINEAR, hllc 

ooo 

ooo 

Cartesian 

400 uni. 

50000KBB 

Qh = 10^9 s-i 

h Jupiter 

RK3, WEN03, hllc 
-i-gravity, -i-thermal cond. 
-i-advection 

fof 

ooo 

Spherical 

240 stretch. 

solar min. 

AttJ - 1315 ergcm^^ s“' 


Notes. Abbreviations: (BC) boundary condition, (f) fixed boundary condition, (o) outflow boundary condition, (SED) spectral energy distribution, 
(uni.) uniform grid spacing, (stretch.) stretched grid spacing, (BB) blackbody spectrum, (solar min.) soI m minimum spectrum ( [Woods & Rottman| 


2002 1 , = Qu/^nr}., Qh - f Lylhvdv, AnJ - f ATrJydv, (RK3) third-order Runge-Kutta scheme (Gottlieb & Shu 


1996 1 , (WEN03) weighted 


essentially non-oscillatory finite difference scheme ^Jiang & Shu||1993} , (hllc) Harten, Lax and vati Leer approxirnate b iemann solver with 
contact discontinuity (|To ro et al.|1994) , (two_shock) nonlinear Riemann solver based on the two-shock approximation (jColella & Woodward] 
[T^IFryxell et al.|2000^. 




Distance (10^® cm) 

Fig. 2. Pure hydrogen weak D-type ionization front. Density (a), tem¬ 
perature (b), velocity (c), and degree of ionization (d) are depicted. A 
CLOUDY simulation (red dashed lines) is compared with a TPCI sim¬ 
ulation (blue solid lines). Initial conditions are given by the gray dotted 
lines. In panel (c) the speed of sound is shown by the black solid line. 
The TPCI simulation agrees well with the independent CLOUDY sim¬ 
ulation. Small oscillations in the high-density region that can be seen in 
panel (c) do not affect the solution significantly. 


photoionization solver, thus, any differences in the results can 
only be caused by the interaction of the hydrodynamics and the 
photoionization solver through the interface. 

Figure shows the simulation of a weak D-type ionization 
front in a pure hydrogen gas in the rest frame of the front to- 
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Fig. 3. Convergence of the models of the weak D-type ionization front. 
The CLOUDY simulation is shown by red, thin lines, the advective 
TPCI simulation by thick, blue lines. The discretization error (dot¬ 
ted), the convergence error (dashed), and the advection length (solid) 
are plotted against the iteration number of the steady-state solver in 
CLOUDY. The grid resolution of the TPCI simulation is indicated by 
the thin, black line. In both simulations the advection length is reduced 
to within a factor of 2 of the grid resolution in the hydrodynamic simu¬ 
lation part. 


gether with the initial conditions. Cold, neutral gas (400 K) with 
a density of 10® cm“^ is irradiated from the right-hand side by a 
hot O-star with a 50 000 K blackbody spectrum. The gas shows 
a steady-state flow toward the source of ionizing radiation. The 
evaporation process is subsonic, starting with 0.1 kms * and 
reaching 7.1 km s ' behind the ionization front. The temperature 
increases to 12 000 K and the density drops to 2 x lO"' cm“^. 

We used TPCI without the advection of species for the ini¬ 
tial settling of the simulation and subsequently restarted the sim¬ 
ulation including the advection as described in Sect. 2.6 The 
settling phase of the dynamical simulation took 3 x 10^* a, after 
which oscillations of the ionization front were reduced to about 
1% of the depth of the ionized region. These oscillations in the 
denser gas are caused by the reflection of waves at the left-hand 
boundary and at the ionization front. The depth of the ionized 
region in the steady-state CLOUDY simulation agrees with the 
dynamic TPCI simulation within the errors produced by the re- 
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Fig. 4. Optically thin cooling function of PLUTO (blue, dotted) and 
CLOUDY (red, solid). The crosses depict the actual cooling in the coro¬ 
nal shock simulation with TPCI. These are a factor of 1.5 lower because 
the simulation shows significant optical-depth effects. 


maining oscillations (see Fig.|^(d)). The width of the ionization 
front is reduced by 25% in the dynamic simulation, which does 
not affect the overall structure, however. 

The convergence of the steady-state solver in the advective 
TPCI and the CLOUDY simulations is shown in Fig.[^ The ad- 
vection length has been iteratively refined whenever the solution 
converged for the current length (see Sect. |2.6| l. Every refinement 
temporarily increases the convergence error, which is then re¬ 
duced over the following iterations until further refinement. The 
final discretization error is on the order of 10%, computed on the 
fine grid of the photoionization solver. This is sufficient because 
the advection length is similar to the coarser resolution of the 
hydrodynamical simulation part. Higher refinement is impossi¬ 
ble because in this case the steady-state solver does not proceed 
beyond this level (see Sect.|^for a discussion of this problem). 

The comparison of the two simulations demonstrates that 
the interface works correctly and is capable of simulating pho- 
toevaporative flows. While the independent CLOUDY simula¬ 
tion with the internal steady-state solver produces smoother pro¬ 
files, TPCI allows studying the dynamical evolution of the phe¬ 
nomenon. 

3.2. Coronal flare - chromospheric evaporation wave 

TPCI can also make use of CLOUDY to compute the cooling of 
collisional plasmas instead of strongly irradiated plasmas. This 
provides the possibility of testing the interface in the optically 
thin limit, where radiative losses only depend on the electron 
number density Ue and the cooling function Ar (e.g., [Sutherland 
|& Dopita|1993] l: 


U^nlAR(T,A). (17) 

The cooling function depends on temperature T and metal abun¬ 
dance A. A radiative transfer solver is not necessary for the 
simulation as long as the optical depth is negligible. Cooling 
in this regime is dominated by emission from ionized metals, 
in which case CLOUDY mostly uses transitions from the CHI¬ 
ANTI database ( |Dere et al.|1997[[Landi et al.|2012| l. The result¬ 
ing cooling function with solar abundances in CLOUDY is sim¬ 
ilar to standard optically thin cooling functions; see Fig.|4 for a 
comparison with the tabulated cooling function that is included 
in the PLUTO code. 

We can thus compare a TPCI simulation with the radiative 
transfer and microphysics solver with a PLUTO simulation using 
the included cooling function. In contrast to the ionization front 


m 

o 

o 

O 




Distance (10 cm) 


Fig. 5. Hydrodynamic evaporation wave in a flaring coronal loop. The 
red dashed lines show a simulation with PLUTO assuming optically 
thin cooling. The blue lines depict the results of TPCI, and the initial 
conditions are given by the gray dotted lines. The density is shown in 
panel (a), the temperature in (b), the velocity in (c), and the cooling 
of TPCI is plotted in panel (d) together with the constant heating rate. 
Panel (d) shows the fractional abundance of the two iron ions with the 
strongest individual contribution to the total cooling rate. In panel (b) 
the temperature decrement in the after-shock evolution is indicated. The 
sound speed is shown by the solid black line in panel (c). Both simula¬ 
tions show the same behavior, except that in the PLUTO simulation the 
cooling rate is twice as strong as in the TPCI simulation, which is due 
to optical-depth effects and small differences of the electron density. 


simulation in the previous section, we now used the same hydro¬ 
dynamics solver (i.e., PLUTO), but different solvers for the ra¬ 
diative cooling. Differences in the results of the test simulations 
can only result from the implementation of the radiative cooling 
from CLOUDY in PLUTO. In this simulation we included all 
30 elements available in CLOUDY, which additionally tests the 
implementation of metals in the inteface. 

We used a setup that was adapted to a strong chromospheric 
evaporation wave in a ID magnetically confined coronal loop. 
Such evaporation fronts appear when a magnetic reconnection 
event releases a large amount of energy in the solar atmosphere, 
which leads to a solar flare (e.g.. Reale & Orlando|2008| l. Evap¬ 
oration of chromospheric plasma into the corona causes a wave 
that propagates at speeds of ~ 400 km s ' through a coronal 
loop, increasing the density by up to one order of magnitude. 
Since optically thin cooling depends on the density squared, the 
cooling rate is increased and, in the long term, the evaporated 
plasma cools down and falls back onto the chromosphere. The 
magnetic loop structure usually survives the flare event. 

Eor the simulation we went into the quasi-rest frame of the 
propagating evaporation wave (see Eig. [^ for the initial condi- 
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tions). The wave front is stable at the start of the simulation, but 
the strength of the wave weakens as a result of cooling. This re¬ 
duces the propagation speed so that during the simulation the 
wave front starts trailing toward the left-hand boundary. This 
proceeds faster in the PLUTO simulation because of the stronger 
cooling rate. Therefore, we do not compare contemporaneous 
states in both simulations but states with the same displacement 
from the initial shock position. There is no irradiation in this 
simulation. 


In the test, coronal material flows in from the right-hand side 
of the domain with a temperature of 3 x 10^ K, a hydrogen num¬ 
ber density of 2.4x 10® cm“^, and a velocity of -780kms“'. The 
plasma is compressed to a density of 7.1 x 10® cm“^ with a tem¬ 
perature of 10® K and a velocity of -300 kms The domain 
length was chosen deliberately longer than a typical magnetic 
loop to emphasize the cooling of the dense gas behind the shock. 

To stabilize a coronal loop model, a constant mechanical 
coronal heating rate has to be applied ( |Rosner et al. 1978 1 . 
This heating of unspecified source balances radiative and con¬ 
ductive losses throughout the loop. We applied a heating rate 
of 5.1/9.8 X 10“^ ergcm^^ s ' in the TPCI/PLUTO simulations, 
which exactly equals the cooling rates in the pre-shock region 
(see Fig. (d)). Advection of species can affect the ionization 
stage of metals behind the shock, which in turn has a weak in¬ 
fluence (factor < 1.5) on the value of the cooling function. The 
increased density affects the cooling rate with a factor of ~ 3^, 
however, which is why we neglected advection of species in this 
simulation. 


Both simulations show the same behavior (see Fig. [^. The 
plasma is heated in the wave front and slowly cools down in 
the post-shock region. The cooling efficiency is a factor of two 
higher in the PLUTO simulation than with TPCI. This has two 
reasons: first, the actual cooling rate in the TPCI simulation is 
lower by a factor of 1.5 than in the optically thin case because 
CLOUDY includes optical-depth effects along the coronal loop 
(see Fig. [^. For example, one of the strongest emission lines 
of the hot plasma is a Fexx line at 12.9 A, which has an opti¬ 
cal depth of 3.9 and does not escape freely. Fe^*® and Fe^** are 
the most abundant iron ions in the hot gas, while Fe^*® and Fe^*^ 
dominate in the cooler gas. Second, the PLUTO simulation com¬ 
putes the electron density on a simple “hydrogen-only” approxi¬ 
mation, which in this case overestimates the electron density by a 
factor of 1.2. Both factors combined result in the total difference 
of the cooling rates (1.2® x 1.5 = 2.2). 

A comparison of the two simulations shows that TPCI can 
also be used for collisional plasmas, especially if optical-depth 
effects are relevant. One drawback is the increase in computa¬ 
tional effort by a factor of 10^ compared to the much simpler, 
optically thin PLUTO simulation. The strong increase in compu¬ 
tational effort results from solving the equilibrium state includ¬ 
ing all metals (see SectQ. For future applications the gain in 
accuracy by using the radiative transfer and microphysics solver 
must be weighed against the increased computational effort. 


3.3. Formation of a Stromgren sphere 

In a third setup, we used the interface to follow the heating of 
cool interstellar hydrogen gas after the ignition of a strong cen¬ 
tral source. This test problem demonstrates the limits of TPCI 
that are due to the underlying assumption of ionization and 
chemical equilibrium. The sudden impact of ionizing photons 
after ignition of the source heats the surrounding gas, and a su¬ 
personic ionization front develops (R-type). The propagation of 
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Fig. 6. Expansion of an R-type ionization front. The TPCI simulation 
(blue solid) is compared with the theoretical evolution of Eq.[T^ Three 
Stromgren radii for different temperatures are indicated by labeled thin 
black lines. The expansion of the ionized region proceeds faster in the 
TPCI simulation because the equilibrium solutions of the photoioniza¬ 
tion solver at each integration step are not restricted by the time since 
ignition of the central source. 


the ionization front is only restricted by the rate of new ioniz¬ 
ing photons emitted by the central source. The rapid expansion 
comes to a rest after recombinations in the ionized sphere bal¬ 
ance the production rate of ionizing photons; the resulting sphere 
is called a Stromgren sphere, and the radius is given by (e.g., |Os-| 
[terbrock & Ferland|2006| l: 


Rs = 


36h 

47TnlaB 


1/3 


(18) 


Here, Qn - Lvjhvdv is the number of ionizing photons emit¬ 
ted by the central source per second, and Ly is the stellar lumi¬ 
nosity per unit frequency interval, no is the initial density and 
ccb = 2.59 X 10 *®cm®s * is the hydrogen B recombination rate 
for a 10000 K plasma. This rate includes recombination to all 
levels but the ground state because recombination to the ground 
state produces an ionizing photon that does not escape locally. 
The R-type ionization front does not raise significant bulk mo¬ 
tion because it proceeds too fast through the interstellar medium 
to cause significant acceleration. 

We started with an initial constant density of 10"' cm~®, 
a temperature of 100 K, and zero velocity. The emission rate 
of ionizing photons is 2 h = 10"'® s ', which is typical for a 
50000 K hot O-star. For the CLOUDY calls we used a spher¬ 
ical grid with the source in the center because geometric dilu¬ 
tion of the radiation field must be included. Since no significant 
bulk flows develop throughout the simulation, we simply used a 
Cartesian grid in PLUTO (irradiation from the center in a spher¬ 
ical grid has not been implemented in TPCI yet), and we also did 
not need to include the advection of species. We only followed 
the time-dependent heating of the initially cool interstellar gas. 

Numerically, R-type fronts have been studied for example by 
Dale et al. ( 2007| l. The authors derived a simple phenomenologi¬ 
cal solution for the expansion of the ionized region after ignition 
of the central source, which can be used to compare the results 
of our code: 


R{f) = Rs (1 - exp {-riQUBt)) 


1/3 


(19) 


This formula is derived by noting that the radius increment per 
time is restricted by the production rate of ionizing photons mi¬ 
nus the recombination rate inside the already ionized region. 

Figure l^shows the propagation of the R-type ionization front 
in TPCI compared with the solution in Eq. 19 Especially the be¬ 
ginning of the simulation shows a faster expansion of the ionized 
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Fig. 7. Temperature evolution during expansion of the R-type ionization 
front. The temperature versus radius is plotted at various times. The time 
step between consecutive lines is constant, the time is referenced by the 
color scheme. The temperature increases from 1000 K to 17 000 K, be¬ 
ing nearly constant throughout the ionized region. The expansion radius 
of the ionized region at a given time is marked by the steep temperature 
decline and increases from 0.085 to 0.16 pc. The distance between con¬ 
secutive lines decreases, indicating the approach to the final Strdmgren 
radius. 



Fig. 8. Thermal conductivity coefficient. The coefficient for electron 
conductivity is depicted by the red dashed line, for neutral hydrogen 
by the green dotted line, and the blue solid line is the sum. The black 
dash-dotted line shows the more accurate weighted sum (Eq. |26^ de¬ 
rived with the actual densities of neutral hydrogen and electrons in a 
CLOUDY simulation without irradiation. In the temperature region be¬ 
low 10 000 K the simple sum is accurate to within a factor of two. 


region than the theoretical prediction. To understand the differ¬ 
ences, we have to remember that CLOUDY seeks equilibrium 
states. Although the evolution seems to be similar to the the¬ 
oretical model, the simulation is based on a different physical 
mechanism. 

Physically, the number of ionized hydrogen atoms cannot ex¬ 
ceed the number of emitted ionizing photons since the central 
source has been switched on. This is not the case in the TPCI 
simulation because the equilibrium photoionization solver as¬ 
sumes that the source already existed for an infinite amount of 
time. The reason that we can nevertheless follow an expansion 
of the ionized region is based on the fact that recombinations 
in the cold gas are more frequent than in an almost fully ion¬ 
ized 10000 K hot region. Therefore, CLOUDY actually com¬ 
putes an equilibrium Strdmgren sphere at every time step, but 
with the temperature structure that is passed from PLUTO (see 
Fig.0. The Strdmgren sphere increases because recombination 
becomes less efficient with increasing temperature. In Fig.j^the 
Strdmgren radii for 5000 K, 10000 K, and 20000 K are dis¬ 
played. At the times when the TPCI simulation reaches the first 
two indicated Strdmgren radii, the mean temperatures in the ion¬ 
ized region are 5600 K and 11 400 K (compare Figs. |^and|7]i. 
The photoionization solver does not assume the simple B recom¬ 
bination rate as in Eq. 18 but allows photons to escape based 
on the optical depth. Hence, recombination in the simulation is 
more efficient, which causes the mean temperatures to be slightly 
higher at the given Strdmgren radii. 

This simulation clearly shows that the user must be aware 
of the restrictions of TPCI and what is actually simulated. The 
interface cannot be used to correctly simulate the expansion of 
an ionized region after ignition of a central source because this is 
an effect of a time-dependent radiation field. It could, however, 
be used to follow the subsequent expansion of the hot gas in 
the Strdmgren sphere into the surrounding cool interstellar gas 
under unchanged irradiation by the O-star because this can be 
approximated by equilibrium states. 


3.4. Simulation of the atmosphere of HD 209458 b 

The intention of the design of TPCI was simulating hot-Jupiter 
atmospheres, and here we present a first test case by solving a 
pure hydrogen atmosphere of HD 209458 b as our fourth test 


case. The expanded atmosphere around HD 209458 b was dis¬ 
covered by |Vidal-Madjar et aT] (|2003|l, and several studies in¬ 


vestigated the system numerically (e.g., |Yelle 2004| Tian et al. 

2005 

Garcia Munoz|2007| 

Penz et al.|2008||Murray-Clay et al. 

2009 

Koskinen et aL|2013 

. These studies provide a solid basis 


fer comparing the results from TPCI. 


For the simulation, we made use of two more modules of 
PLUTO: gravitational acceleration and thermal conductivity. 


3.4.1. Gravity 

The simulation of the escaping hot-Jupiter atmosphere incorpo¬ 
rates the gravitational and centrifugal acceleration in the two- 
body system, also referred to as tidal forces. With the planet in 
the center, we derive for the effective potential 


(l)(r) = 


GMpi 


a — r 


■ 2/, 

-CO (/cm ■ 


rf. 


( 20 ) 


Here r denotes the radial distance on the axis toward the star, 
Mpi and M^t are the planetary and stellar mass, a is the semi¬ 
major axis, 1cm is the distance to the center of mass, and G is 
the gravitational constant. The angular frequency w is given by 
Kepler’s third law 


(Inf 2 G(Mst + Mpi) 


( 21 ) 


In PLUTO we include the acceleration: 



GMpi GMst 

(a - r)2 


G(Mst + Mpi) 


(^cm r'). 


( 22 ) 


By plotting the equation with the values of the system 
HD 209458 ( [Wright et al.|2011] ), we find that on the star-planet 
axis material can be bound to the planet up to 4.18 Rpi, which is 
the size of the planet’s Roche lobe on this axis. 


3.4.2. Thermal concductivity 

We used a simplified approach to include thermal conductivity 
by adding the conductivity of an electron proton plasma and of 
a neutral hydrogen gas. The first dominates for temperatures in 
excess of 50 000 K, the latter below this. 
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The following coefficients are taken from |Banks & Kockarts| 
( 1973| l. The heat flux is given by 

Fc^-ttVCT). (23) 


For highly ionized gases, thermal conductivity is dominated by 
electron-ion collisions ( Spitzerfl 978) 1. The conductivity coeffi¬ 
cient is given by 




(24) 


The thermal conductivity coefficient of neutral hydrogen is pa¬ 
rameterized as 


/thh = 379 ® 


(25) 


and the total coefficient is given by the weighted sum 

/<• = - V = - (Wetfe + nuKii) . (26) 

n n 

In our simulations we used the direct sum of the two coefficients: 


■kh- 


(27) 


In the relevant temperature range below 10000 K the approxi¬ 
mation is correct to within a factor of two (see Fig. |^. Since we 
do not expect thermal conductivity to have a significant impact 
on the results ( Garcia Munoz||2007 1, the procedure is sufficient 
to verify this claim in our simulation. 


3.4.3. Simulation 

We used a pure hydrogen atmosphere, irradiated from the top 
by the host star, and simulated the escape at the substellar point. 
With a stretched grid the cell size is increased from 0.0005 Rpi 
in the lower atmosphere to 0.087 Rpi at the top. The temperature 
in the lower atmospheres was fixed at the planet facing boundary 
to 1000 K, which is approximately the average dayside bright¬ 
ness temperature of HD 209458 b ( jCrossfield et al.|2012| l. In the 
simulation this was implemented by fixing the lower boundary 
density to lO'^ cm~^ and the pressure to 140 dyn cm“^ , which is 
equivalent to 0.14 mbar. The gas is heated by ionizing radiation 
with a solar minimum EUV spectral energy distribution (Woods 
& Rottman|2002|l, which is a good approximation for an inactive 


GO-type dwarf ( |Sanz-Forcada et al.|2011[|Montes et al.|2001| l. In 
the spectral range below 912 A (XUV), the planetary atmosphere 
at the distance of 0.047 AU from the host star is irradiated with 
a flux of Fxuv = 1315 ergcm^ 
to previous estimates, for instance, 1800 erg cm 


'. The derived value is similar 
' iKoski- 


nen_et al.1(|2013|l. For this inactive host star (|Sanz-Forcada et aTT 


2011 1, X-rays (< 100 A) contribute only 3% of the energy in the 


XUV spectral range. In contrast to active host stars. X-rays have 
no strong impact on the mass loss of HD 209458 b. The initial 
conditions of the simulation are shown together with the results 
in Fig.[9| A comparable setup can be found in the simulations of 
|Penz et al.|(|2008|l, |Murray-Clay et al.|(|20091l or|Koskinen et al.| 
02013^ . 

In the expanding atmosphere neutral hydrogen is transported 
from the lower to the upper atmosphere, which means that advec- 
tion of species is essential for the results in this simulation. To 
save computational time, we ran the simulation without the ad- 
vection of species until a steady-state planetary wind was found 
and all shock waves had subsided. The resulting atmosphere was 
used as initial atmosphere for the final simulation including the 
advection of species. The convergence of the advective scheme 
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Fig. 9. Evaporating hydrogen atmosphere of HD 209458 b. Density (a), 
temperature (b), velocity (c), ionization fraction (d), heating and cool¬ 
ing rate (e), and the heating fraction (f) are plotted versus the planetary 
radius. The red dashed lines depict the initial settling phase, the blue 
solid lines the final TPCI simulation including advection of species. 
The initial conditions are given by the gray dotted lines. Panel (a) also 
shows the density structure of an isothermal atmosphere with 500 K 
(black solid line). In panel (c) the speed of sound is indicated (black 
solid line), and in panel (e) the three main cooling agents in the atmo¬ 
sphere are shown by black thin lines. “H line emis.” is cooling due to 
line emission from hydrogen atoms (mainly Lya), “ff emis.” is cooling 
by free-free emission of electrons, and “H recomb.” is recombination 
cooling of hydrogen. Additionally, thermal conductive cooling (cond.) 
is shown by a black dotted line. The atmosphere expands in a transonic 
flow. A small amount of the absorbed energy (~ 11%) is lost through 
free-free and Lyo- emission. 


does not proceed beyond an advection length of 0.27 Rpi. This is 
comparable to the grid resolution at the top of the atmosphere, 
but much higher than the resolution in the lower atmosphere, so 
that effects due to the advection of species are initially only re¬ 
solved in the upper atmosphere. However, below 1.05 Rpi veloci¬ 
ties are low and the advective timescale is longer than the hydro- 
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gen recombination timescale. Advection of species can safely be 
neglected below this level. Furthermore, we manually reduced 
the advection length to 0.01 Rpi and did not find significant 
changes in the solution. We therefore conclude that the advec¬ 
tion of neutral hydrogen into the upper atmosphere is resolved 
throughout this simulation. 


3.4.4. Results 

The atmosphere expands in a steady-state transonic flow that is 
directed from the hot Jupiter toward the host star. We find a steep 
density and temperature gradient within the hrst 0.2 Rpi (see 
Fig.0. At 1.08 Rpi the atmosphere shows a temperature mini¬ 
mum with only 500 K at a pressure of 0.7 //bar; the density strat- 
ihcation below this point differs only slightly from an isothermal 
atmosphere with the same temperature (see Fig.[^(a)). EUV ra¬ 
diation does not penetrate to this depth, but is mainly absorbed 
from 1.1 to 1.2 Rpi, where the temperature of the atmosphere 
strongly increases. In this pure hydrogen atmosphere free-free 
emission is the main radiative cooling agent below 1.08 Rpi, 
reemitting about 50% of the total heat input, up to 70% of which 
is provided by hydrogen line absorption at this depth. Hydro¬ 
gen line emission, mainly Lya radiation, dominates the radiative 
cooling from 1.08 to 2.0 Rpi, and recombination of hydrogen is 
the most efficient radiative cooling agent at higher altitudes. 

Averaged over the atmosphere above the temperature mini¬ 
mum, 89% of the heat input is used for heating and acceleration 
of atmospheric gas, and eventually for lifting the material out of 
the gravitational potential well of the planet. The total energy 
gain is divided as follows: 77% of the radiative energy input is 
converted to gravitational potential energy, 17% to kinetic and 
6% to internal energy. Note that we dehne the heating fraction 
as 


/h - 



(28) 


here Gr is the actual heating produced by the absorption of ra¬ 
diation, not the energy of the absorbed radiation (see Eq. |^. In 
our hydrogen-only simulation, only free-free and Lya emission 
decrease the radiative heating fraction signihcantly (see Eig. 
(f)). The fraction approaches 1.0 in the upper atmosphere; how¬ 
ever, emission by helium or metals is not included in this sim¬ 
ulation. Especially line absorption and emission by metals can 
significantly affect the heating and cooling rates in the lower at¬ 
mosphere. 

Replacing the denominator in Eq|^with the absorbed radia¬ 
tive energy gives the net heating efficiency, which is often used 
in hydrodynamic simulations of hot-Jupiter atmospheres to ac¬ 
count for the fraction of energy that goes into ionizing hydrogen 
and also for not specihcally included cooling agents. In this par¬ 
ticular simulation, stellar emission is absorbed from 50 to 912 A 
by the atmosphere (more than 50% absorption). This absorption 
produces 92% of the total heating rate; the remaining 8% are 
caused by hydrogen line absorption, mainly Lya absorption. The 
heating efficiency in the XUV range is 0.71. 

Thermal conduction does not play a significant role in the at¬ 
mosphere since it is more than two orders of magnitude smaller 
than the main cooling agents at every point throughout the at¬ 
mosphere (see Eig.[9](e)). Comparing the structure from the ini¬ 
tial settling phase without advection of species with our hnal 
results shows that the advection of neutral hydrogen into the up¬ 
per atmosphere increases the average temperature from 3000 to 
5500 K. Advecting neutrals into the otherwise highly ionized 


region increases the heating rate because more ionization pro¬ 
cesses occur. This also boosts the heating fraction, as can be seen 
in Eig.[^ The final velocity exceeds the sonic speed of 10 km ’. 

Our results are in excellent agreement with the C3 model 
of Koskinen et al.| ( T013| l for the atmospheric evaporation at the 
substellar point. Assuming a spherically symmetric atmosphere, 
we derive a mass-loss rate of M - Anr^pv - 1.4 x lO" gs“', 
which is equal to 0.0023 Mjup Ga *. This value is an upper limit 
since we simulated the substellar point with the highest mass- 
loss rate. Eor the surface averaged mass loss, a factor of 1/4 was 
used by |Koskinen et al. which places our result just slightly be¬ 
low their values, ranging from 4 - 6.5 x 10'° gs“'. Within 5% 
accuracy, advection of species does not change the mass-loss rate 
because a slightly higher velocity in the advective simulation is 
canceled by a lower density. 

|Murray-Clay et al/] ( |2009| l were the hrst authors to include a 
model for Lya cooling in the simulation of escaping planetary at¬ 
mospheres. Their spatial distribution and total cooling rate com¬ 
pares well with our results. Lya cooling leads to a lower peak 
temperature in the atmosphere and reduces the mass-loss rate 
( [Koskinen et al.|2013| l. Note, however, that deeper in the atmo¬ 
sphere the density is high enough that collisional de-excitation 
dominates and absorption of Lya radiation is a heating agent. 
About one quarter of the total cooling through Lya radiation is 
counterbalanced by Lya heating in the lower atmosphere. On the 
other hand, Murray-Clay et al. used a gray absorption scheme, 
which leads to a mass-loss rate lower by a factor of four ( jKosk-j 
jinen et al.|2013| l. Such differences clearly prove the advantages 
of the detailed photoionization solver included in TPCI. Neglect¬ 
ing any significant cooling or heating agent affects the mass-loss 
rate, which crucially depends on the available energy. This be¬ 
comes even more important when metals are included in the at¬ 
mosphere and an increasing number of processes must be mod¬ 
eled. 

Although the mass-loss rate does not change when the ad¬ 
vection of species is included in the simulation, the advection of 
neutral hydrogen into the upper atmosphere is important for an¬ 
alyzing observational results. The Lya absorption signal of the 
planet during transit does change signihcantly because of the 
higher neutral hydrogen column density. The high velocities of 
the expanding atmosphere are crucial to transport the neutral hy¬ 
drogen into the upper atmosphere. 


4. Discussion and conciusion 

We have presented TPCI, which is an interface between the 
MHD code PLUTO and CLOUDY, an equilibrium photoioniza¬ 
tion solver. TPCI enables simulations of hydrodynamic steady- 
state and slowly evolving hows, including multiple physical ef¬ 
fects such as gravity and thermal conduction. It simultaneously 
solves the equilibrium state of a gas or plasma under strong ir¬ 
radiation. The code is valid from cold molecular regions to fully 
ionized hot plasmas and computes the radiative transfer along 
ID rays. 

The presented test simulations prove the validity of the com¬ 
bined simulations. In the introduction we have established a 
number of requirements for TPCI. The interface has proven its 
capabilities in solving the ionization and microphysical state of 
gases for a wide range of densities and temperatures. The second 
test problem (see Sect. |3.2| i demonstrates the simulation of gases 
including metals and shows effects of X-ray radiative transfer. 
Multidimensional simulations or magnetic helds have not been 
studied in our verihcation of the ID scheme. In general, we hnd 
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Fig. 10. Execution times of CLOUDY runs on local cpu with setups 
based on the hot-Jupiter simulation. Simulation times vary in the shaded 
range based on density and temperature. Each step increases the com¬ 
plexity of the simulation, starting with a simulation including only hy¬ 
drogen and 10 zones in CLOUDY. In this case the execution time is 
dominated by the initialization of CLOUDY. The next steps are with 
100 and 1000 zones. We then include step by step molecules, helium, 
carbon, nitrogen, and oxygen, and all 30 lightest elements. Einally, we 
iterate twice in CLOUDY over the complete domain to correct optical- 
depth effects. 


that TPCI runs stably and reliably. The combination of the adapt¬ 
able hydrodynamics code with a versatile radiative transfer and 
microphysics solver is generally applicable to photoevaporative 
flows and can for example be used to study the dynamics of H it 
regions or the evaporation of protoplanetary disks. Furthermore, 
CLOUDY is a powerful tool for spectral analysis, and observed 
spectra can be directly compared with results from the simula¬ 
tions. 

One drawback in the use of TPCI is the computational effort 
of the detailed photoionization solver, which usually consumes 
almost all of the computational time. For the simulation of the 
test problems we used a few days to up to two months on a single 
3.1 GlTz processor. One-dimensional simulations are serial, so 
that simulations do not perform faster on clusters. In the search 
for hydrodynamic steady-state situations, the best call for reduc¬ 
ing the computational effort is to hnd good initial conditions to 
speed up the convergence and to include more complex physics 
step by step from pre-converged simpler solutions. 

Figure [1^ shows execution times in single calls to the pho¬ 
toionization solver. The setup for the calls is based on the simu¬ 
lation of the hot-Jupiter atmosphere, but we used constant hy¬ 
drogen number densities ranging from 10^ to lO'^ cm~^ and 
constant temperatures ranging from 1000 to 10000 K. The ac¬ 
tual execution times vary depending on the parameters. In a 
typical simulation of a hot-Jupiter atmosphere CLOUDY needs 
about 1000 zones for the adaptively adjusted grid, which means 
that execution times were about 100 s, but two iterations are 
needed to correctly compute the optical depth for the escape 
probability formalism. The settling phase of the simulation took 
about 10000 calls to the photoionization solver, which explains 
the simulation time of about one month. For the advection of 
species, at least 30 iterations have to be made in the CLOUDY 
internal steady-state solver, but when starting from the previ¬ 
ously converged simulation, about 200 calls to the photoioniza¬ 
tion solver suffice, which adds up to one day. In the future, the 
abundance of metals can be linearly increased during the simu¬ 
lation to prevent strong oscillation in the atmosphere and, thus, 
limit the CLOUDY calls to a minimum. 

The convergence of the advection of species in the CLOUDY 
internal steady-state solver can cause problems in certain situ¬ 
ations. The result is a solution where effects of the advection 
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of species are only resolved on scales larger than the advection 
length. One possible solution is to manually reduce the advec¬ 
tion length to the needed accuracy, which increases the compu¬ 
tational effort, however, because more iterations in CLOUDY 
are needed. Furthermore, advection can only be included in ID 
simulations with flows toward the irradiating source. A possi¬ 
ble solution would be to include passive scalar helds in PLUTO 
for every species and solve the advection by explicitly including 
source and sink terms, which could be retrieved from each call to 
the photoionization solver. This complicates the interface, how¬ 
ever and necessitates more essential changes of the CLOUDY 
source code because the number densities of each species would 
have to be passed to CLOUDY and source and sink terms must 
be retrieved. Species include all ionization stages of the 30 light¬ 
est elements plus all included molecules. 

The simulation of atmospheric evaporation in the well- 
studied hot Jupiter HD 209458 b demonstrates the capabilities of 
TPCI. The simulation of the hydrogen-only atmosphere shows 
significant differences compared to |Murray-Clay et aL] ( |2009| l 
because the authors used a gray absorption scheme, or com¬ 
pared to Penz et al. ( 2008| l, who neglected Lya cooling. Our at¬ 
mosphere shows an excellent agreement with the corresponding 
model of Koskinen et al. ( 2013p . We resolved a temperature min¬ 
imum in the lower atmosphere because our simulation continues 
deeper into the atmosphere. Furthermore, we used a more com¬ 
prehensive microphysics solver, which reveals significant cool¬ 
ing due to free-free emission in the lower atmosphere, which has 
previously been underestimated ( |Murray-Clay et al.|2009[ ). 

The increased accuracy in solving the absorption and emis¬ 
sion processes is a clear advantage, and by using the well-tested 
CLOUDY code, we did not neglect any signihcant processes in 
our model. The total change in the mass-loss rate compared to 
the model of [Koskinen et al.| ( |2013| l is lower than a factor of two 
in this pure hydrogen atmosphere, but larger differences can be 
expected when metals will be included. In the future, we will 
extend this initial test case to include molecules and metals and 
will investigate the effects of strong X-ray irradiation in the at¬ 
mospheres of planets around active host stars. 
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